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Abstract 

Recursive maps of high order of convergence m (say m = 2 10 or 
m = 2 20 ) induce certain monotone step functions from which one can 
filter relevant information needed to globally separate and compute the 
real roots of a function on a given interval [a, b]. The process is here called 
a root distiller. A suitable root distiller has a powerful preconditioning ef¬ 
fect enabling the computation, on the whole interval, of accurate roots of 
an high degree polynomial. Taking as model high-degree inexact Cheby- 
shev polynomials and using the Mathematica system, worked numerical 
examples are given detailing our distiller algorithm. 


1 Introduction 

By a higher-order root distiller we mean an algorithm to compute, simultane¬ 
ously, all (or almost all) the real roots of a polynomial / of (high) degree d, in a 
given interval [o, b]. The algorithm relies on a single application, on the interval, 
of a map g of (very) high order of convergence m (for instance m = 2 10 = 1024 
or m = 2 20 = 1 048 576). The approach may be seen as a modern computational 
perspective of the global Lagrange’s ideal pQ. 

Once defined such a higher-order map g, the roots of the function / in the 
interval are obtained from a table C = {(x\,yi), (x 2 , IJ2), ■ ■ ■ , (%n, Vn )}> where 
Hi = g{xi), and Xi belong to an uniform grid (of IV+ 1 nodes) of width h, defined 
on [a, b\. For suitable choices of h and k (the parameter k controls the order of 
g ), the map g induces an invariant monotone step function leading to certain 
subsets of C, say S\, S 2 , Si,S r . These subsets will be called ‘platforms’ 
(see Section |3| and have the property that the second component of the points 
(x , y)i on the platform Si are equal to (or close to) the i-th root of the polyno¬ 
mial. 

A graphical inspection of the list C may be useful not only to observe the 
distribution of the roots in [a, b\ but also to suggest good choices for the two 
parameters controlling the algorithm, which are the mesh width h and the index 
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k, the latter related to the order of convergence m of the map g (m = 2 k+1 in 
the case of simple roots). 

Although in this work we only deal with roots of polynomials, the same proce¬ 
dure can be adapted to non-algebraic equations f(x) = 0 having at least one 
root in a given interval, or with the case of multiple roots or even to functions 
in R n 0. 

In finite arithmetic, one of the main feature of our distiller process is that, by 
construction, the values yi are generally quite immune to rounding error prop¬ 
agation. Therefore a roots’s distiller can be seen as a powerful pre-conditioning 
instrument, in particular for polynomials whose coefficients are numeric. The re¬ 
ferred immunity to rounding error propagation is closely related to the fact that 
a map of higher-order of convergence leads necessarily to a stationary ‘monotone 
machine step function’ — if the recursive process which generates the map g is 
taken appropriately, that is, for k sufficiently large. Details on the monotone 
machine step functions are further explained in sections [2] and [3] 

Thanks to the super-attracting property of a map of a high order of convergence, 
each tread (or ‘platform’) of the monotone machine step function - correspond¬ 
ing to the theoretical subsets S. t referred above- contains several machine accu¬ 
rate values yi which are close approximations of the zeros of the given function 
/. In general, all the necessary information in order to approximate the zeros 
of a given map with a prescribed accuracy is contained in these treads. 

Our root distiller is constructed in order to overcome some common numerical 
issues appearing in the computation of zeros of a given function and in partic¬ 
ular of roots of polynomials of high degree. It is well known the inherent ill 
conditioning of the computation of polynomial roots. For instance Mathematica 
commands for approximating roots of polynomials of high degree may produce 
useless numerical results when low-precision finite arithmetic is used. On the 
other hand, dealing with exact polynomials of high degree d, say d > 100, pre¬ 
vents us from using exact arithmetic due to CPU excessive cost. 

To be more precise, suppose that a numeric expression for the (first kind) Cheby- 
shev polynomial of degree 40 is defined by the command N [ChebyshevT [40, 
x ], 8], where the coefficients are deliberately forced to have 8-digits precision. 
The commands Solve, Reduce, and Roots produce useless numerical results (cf. 
paragraph 1 1.1[ ) since the computed roots are heavily contaminated by rounding 
error (even though the degree d = 40 of such polynomial is moderate). Our root 
distiller deals efficiently not only with this case but it also produces accurate 
answers, for instance with a 500-degree Chebyshev polynomial. 

In Section [2] we detail the construction of a specific map g. For that, it is given 
a positive integer prec and two parameters h and k. The parameter k controls 
the order of the map g to be constructed, h is the mesh size and prec fixes the 
precision to be used in the computations of the images yi = g(xi) in the list C. 

In Section [3] it is illustrated how a map g of high order of convergence leads to a 
monotone step function which contains the relevant information to be distilled. 
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We chose as basic model a 4-degree Chebyshev polynomial of the first kind with 
prec = 8 and the parameters h = 0.1 and k = 3. The respective map g has 
order of convergence 16 and the absolute error of computed roots is of order 
1CT 8 , meaning that the accuracy used on C is preserved. 

Mathematica code is presented in sections[2]an <10 including the process used for 
filtering the relevant values in the respective list C (other filtering possibilities 
may also be considered). 

Numerical examples have shown the efficiency of the proposed distillers. In 
particular, we construct here a distiller for the computation of the positive roots 
of the Chebyshev polynomial of degree 500, defined in [—1,1], with precision 
forced to be prec = 5000 and parameters h = 0.00025 and k = 20 (and so the 
respective map g has order of convergence m = 2 21 = 2 097 152). The computed 
roots have 5000-correct digits. 

An automatic choice of appropriate parameters h and k in order to achieve 
a preassigned tolerance error can be done, but this is out of the scope of the 
present work. 

1.1 Motivation: a low precision Chebyshev polynomial 

Setting the precision prec = 8, we obtain the following Mathematica expression 
for the Chebyshev’s polynomial of degree 40, N[ChebyshevT[A 0, a;], prec]: 


1.0000000 - 800.00000a; 2 + 106400.00a; 4 - 5.6179200 x 10 6 a; 6 + 

+1.5690048 x 10 8 x 8 - 2.6777682 x 10 9 a; 10 + 3.0429184 x 10 lo a; 12 - 
-2.4343347 x 10 n a; 14 + 1.4240858 x 10 12 a; 16 - 6.2548083 x 10 12 a; 18 + 

+2.1002988 x 10 13 a; 20 - 5.45 5 32 1 5 x 10 13 a: 22 + 1.102 9 2 37 x 10 14 a; 24 - 

-1.7375290 x 10 14 a; 26 + 2.123 6 4 66 x 10 14 a; 28 - 1.9918340 x 10 14 a; 30 + 

+1.4055280 x 10 14 a; 32 - 7.2 1 55 451 x 10 13 a; 34 + 2.5426206 x 10 13 a; 36 - 

-5.4975581 x 10 12 a; 38 + 5.4975581 x 10 n a; 40 . 

The Mathematica commands Solve, Reduce and Roots produce, respec¬ 
tively, the following useless output: 

{-0. X 1CT 1 , -0.35, -0.271, -0.1951, -0.11754, -0.0392598, In . x 1 0+ r 

In - x i o -1 l. In . x i n-+ In . x i EZZ1 , In . x i n-+ In . x i Q-+ In . x i IQ, In . x i n-+ In . , i q_3 . 

In . x i n-+ In . x i n-+ In . ■ 1 53 , In . x i n^l. In . x i n-+ In . x i n+ In . x 1 53 , In . ■ 1 53 , 

In . X 1 53 , In . X 1 53, In - x 1 53, In . x i n+ In . x i n 1 ! . In . x i n 1 ! , In . x 1 53 , In . ■ 1 53 , 

In . x i o+ In . ■ i n+ In . x i El, 0 . 0392598 , 0 . 11754 , 0 . 1951 , 0 . 271 , 0 . 35 , o.xio -1 } 

x = -0.4 I I x = -0.35 I | x = -0.271 I I X == -0.1951 | | x = -0.117537 | | x = - 0.0392598 | ] 
x = 0 | | x = 0.0392598 | | x == 0.117537 | | x == 0.1951 [} x = 0.271 | | x = 0.35 | | x == 0.4 


X == 0 I I X == 0 I I x == 0 I I X == 0 I I x == 0 I I X == 0 I I X == 0 I I X == 0 I I X == 0 I I X == 0 I I 

x == 0 | | x == 0 | | x == 0 | | x == 0 | | x == 0 | | x == 0 | | x == 0 | | x == 0 | | x == 0 | | x == 0 | 

x == 0 | | x == 0 | | x == 0 | | x == 0 | | x == 0 | | x == 0 | | x == 0 | | x == 0 | | x == 0.0392598 

x == -0.0392598 | | x == 0.11754 | | x == -0.11754 | | x = 0.1951 | | x == -0.1951 | | 

x == 0.271 | | x == -0.271 | | x == 0.35 | | x == -0.35 | | x == 0. x 10 _1 | | x == - 0. x 10 _1 
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Of course the commands NSolve and NRoots also give useless values. 

We aim to obtain ‘machine’ acceptable answers, that is to compute the real 
roots of the 40-degree Chebyshev polynomial in the interval [a, b\ = [—1,1], 
which are simple and distinct, with an accuracy close to that of the data (recall 
that 8-digits precision has been assigned to the polynomial coefficients). 

The algorithm which we call ‘root distiller’ is described in what follows, and 
shows to be able to accomplish such a desideratum. The code can easily be 
included in a single function in order to produce the referred machine point list 
C. once predefined the function f, the bounds of the interval, the preassigned 
precision prec and the mesh size h. After a convenient filtration of the data 
in CL. the respective output should be considered global in the sense that it is 
able to (simultaneously) produce accurate approximations of the roots in the 
interval, as well as realistic error estimates to each of them (see Section 4). 

2 Higher-order educated maps and monotone 
step functions 

In a global approach to roots’s computation by means of a smooth high order of 
convergence map j:[o,i]d —> R, many of the domain points are irrelevant, 
in the sense that their image under g might not be a number or is repealed 
from a fixed point of g. In fact any map of order of convergence greater than 
one enjoys such a repealing/attracting property - like in the well known cases 
of the Newton’s or secant methods for approximation of simple roots. So, once 
defined a map g of sufficiently high order of convergence, the points which are 
not a number, nor in the interval [a, b ], neither attracted to a fixed point will be 
ignored. This is the reason why we then will call g an ‘educated’ higher-order 
map. 

In general, the recursive maps to be considered have order m of convergence 
which can go up to m = 2 10 = 1024, or greater. The recursive process used to 
define the map makes possible to obtain a monotone step function defined in 
the interval [o, b]. From this step function one extracts the relevant computed 
g -images through a filtering process in order to obtain as output most, or all, 
the roots of the equation f(x) = 0. In particular, our distillers will allow us 
to compute the roots of a Clrebyshev’s polynomial of high degree, a task not 
feasible by the exact methods provided by the Mathematica system (version 
10.02.0 running on a Mac OS X personal computer has been used in this work), 
unless the interval is small and the working precision high. 

2.1 Recursive construction of the higher-order map g 

We now explain the recursive construction of a map g of high order of conver¬ 
gence by taking as a seed the Newton’s map. The map g is the fc-fold composi¬ 
tion of this seed and has order of convergent m = 2 k+1 . The construction of g 
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goes through and ‘education’ process aiming to obtain a map satisfying a fixed 
point theorem in the interval [a,b]. More precisely, g is constructed in order to 
satisfy the following properties: 

(i) g([a,b]) C [a,b] . 

(ii) | g(x) — x\ ^ (b — a), for all x for which NumericQ [g(x)l is True . 

(iii) The points x € [a, 6] not satisfying (i) and (ii) are ignored (a Null is 
assigned to g(x)) . 

For k sufficiently high, the educated map g will act on [a, b] as a kind of a ’mag¬ 
net’ having both good theoretical and computational properties. This ‘mag¬ 
netic’ property is better perceived by inspecting a plot of the respective induced 
monotone step function. 

It can be proved that for a fixed mesh size h, an educated map induces an in¬ 
variant (or stationary) monotone step function, whenever the folding parameter 
k is sufficiently large. This invariant step function will be called the machine 
step function associated to a map g = g[h, k\. 

Choosing suitable values for k, the second component of points on the treads 
or ‘platforms’ of the associated step function contain (by construction) accurate 
approximations of the roots of f(x) = 0. Moreover, the platforms of such step 
function are automatically sorted in increasing order of their heights, defining 
so a monotone step function in the interval [a, 6]. The later filtering process of 
the data of this step function will hopely solve the referred global Lagrangian 
root’s problem. 

In the following illustrative example an uniform mesh of points, of width h , is 
defined on the domain range [a, b\. The ListPlot command is used in order to 
observe the behaviour of an higher-order educated map g on the referred mesh. 

Although in this work the seed used in the recursive process is the Newton’s 
map, any other method of order of convergence greater than one could be used. 
For instance, the secant method, of order (l + v / 5)/2, and Ostrowsky’s methods, 
of orders > 3, are other obvious options. 

2.2 The map g from the Newton’s seed 

Fixing a precision prec, assume that the numeric expression for a given function 
/ is in memory as well as the bounds a and b of the interval where the roots 
of f(x) = 0 are required. Given the (folding) parameter fc, the following code 
defines a general recursive function g[x,prec] using Newton’s map as seed (see 
below functions newt [0,x,prec] and its recursive version newt [k,x,prec]). 
The map g , of order of convergence m = 2 k+1 , is given below as the function 
named g[x,prec]. 

Note that when y=g[x,prec] is a number, the assigned precision to y is forced 
to be the same as the precision of x. This prevents the Mathematica system 
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to correct the output of each calculation of y in the case it occurs of a loss of 
significant digits. The code for the function g[x,prec] follows. 

(* Messages such as infy or indet are ignored: their occurence *) 

(* means that the map g has a singularity at a point x in [a,b] *) 

Off[Power::infy]; 

Off[Infinity::indet]; 

(* The function newt[0,x,prec] *) 

(* gives the standard Newton's image of the point x. Both pre-image *) 

(* and image are forced to have prec digits assigned. *) 

newt[0, x_, prec_] := newt[0, x, prec] = Block[{xx = SetPrecision[x, prec]}, 
SetPrecision [xx - f [xx] / f ' [xx] , prec] ] ; 

(* For the non-negative integer k,the function newt[k,x,prec] gives *) 

(* the k-fold composition of Newton's map. *) 

newt [k_, x_, prec_] :=newt[k, x, prec] =newt[0, newt[k-l, x, prec], prec]; 


(* The higher-order map g[x,prec] gives the image *) 

(* for the k-fold function newt[k,x,prec]. *) 

(* This map is educated in order that the image y lies in [a,b]. *) 

(* To an image not in [a,b] or whose step=y-x is not in [a,b] it is *) 
(* assigned the Null symbol. *) 


g[x_, prec_] := Block[{y, step}, 
y = newt[k, x, prec]; 
step = y - x; 

(* capture educated numeric points in [a,b] : *) 

If [Abs [ step] <= (b-a) && a^y^b , y] 

i ; 


3 An illustration with a low precision 4-degree 
Chebyshev polynomial 

As an illustration of the occurrence of a monotone step function induced by 
the map g , let us consider a 4-degree Chebyshev polynomial in the interval 
[a, b} = [—1,1], and define an uniformly spaced mesh of width h = 0.1 (that is 
21 equally spaced nodes). We apply the above g[x,prec] code, with assigned 
parameters prec = 8 and k = 3, that is, in this case the map g has order of 
convergence m = 2 k+1 = 16. 

In Figure[l]it is displayed the plot of the respective point list C and the computed 
values g[x,prec\. This figure is self-explanatory: there are 4 roots corresponding 
to the 4 platforms in the displayed graphic; an increasing step function is sug¬ 
gested by the dotted broken line. Each of the four observed platforms is formed 
by a set of points and from this set of pairs (x, y) we filter the value y of the pair 
which has the second component closer to the first. The value of y filtered this 
way is an approximation of a polynomial’s root. In this example, the 4 points 
filtered are displayed in the last column of the table. These 4 values are 8-digit 
accurate roots of the Chebyshev polynomial of degree 4. 

Looking at the second column of the table in Figure 1, it is clear the super- 
attracting property of the ‘magnet’ g: all the numeric points in each platform 
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solutions 


18 super-attracted points 



Xi 

g[xi,precj 

-1. 

-0.92387953 

-0.9 

-0.92387953 

-0.8 

-0.92389110 

-0.7 

out 

-0.6 

-0.38268343 

-0.5 

-0.38268343 

-0.4 

-0.38268343 

-0.3 

-0.38268343 

-0.2 

-0.38268343 

-0.1 

0.92388432 

0. 

out 

0.1 

-0.92388432 

0.2 

0.38268343 

0.3 

0.38268343 

0.4 

0.38268343 

0.5 

0.38268343 

0.6 

0.38268343 

0.7 

out 

0.8 

0.92389110 

0.9 

0.92387953 

1. 

0.92387953 


-0.92387953 


-0.38268343 


0.38268343 


0.92387953 


Figure 1: A step function associated to a 8-digits precision, 3-fold map g , applied 
to a 4-degreed Chebyshev polynomial. 


have the second component very close to the respective exact fixed point of 
the map g. Repetting the computations for k = 4, the corresponding table 
is identical to the one in Figure 1, except the image of x = —0.8, which is 
—0.92387953. This means that the respective map g , now of order 2 k+1 = 32, is 
invariant and so the former computed roots have indeed 8 correct digits. Saying 
it in other words - the 4 computed values are ‘machine’ fixed points for this 
map (see the last column in Figure 1). 


4 Distilling Chebyshev polynomials of high de¬ 
gree 

For a given precision prec and appropriate parameters h and k, we developed 
a simple Mathematica code which will be tested in order to approximate the 
roots of high degree Chebyshev polynomials. A realistic estimative of the error of 
each computed root is also easily obtained. In fact, if g has order of convergence 
m > 1 and y = g(x) is a value close to a fixed point a , the error of y satisfies 

« - y - g(y) ~ y- 

Concerning our polynomial models, since the roots become closer when we in¬ 
crease the degree of the Chebyshev’s polynomial, the value of the mesh size h 
and the parameter precision prec need to be adjusted accordingly. Our aim is to 
obtain a highly accurate bound for the positive roots of a 500-degree Chebyshev 
polynomial of the first kind. 


7 


















4.1 Example: a 500-degree Chebyshev polynomial 

A 500-degree Chebyshev’s polynomial is highly oscillating and so its roots are 
very close. Therefore, it is necessary to set a sufficiently large value of the 
precision prec, and choose a convenient mesh size h in order to obtain the 
required numerical results. We consider now the interval to be [0,1]. 

In order to observe the platforms of our distiller, we display the plots of the 
educated maps g for some values of the parameter k and of prec as a guide for 
the choice of the right values of these parameters. The next figure compares the 
graphics of the function g [x, 100], respectively for k = 0 and k = 10 (the mesh 
size in the plot is not uniform since it is automatically generated by function 
Plot in the interval [a, b]). 

It is clear from Figure [2] that a 100-digits precision is not enough to obtain 
the roots localised on the right-half domain, while one can expect a good sep¬ 
aration of roots in the interval by setting prec = 5000, as suggested by Figure[3] 




Figure 2: For prec = 100 and a mesh automatically generated by Plot , only the 
roots localised on the left half-domain are detected. 


An advantage of choosing an higher order of convergence map g will become 
more apparent if we restrict the interval to [0.9,1.0], which contains the desired 
greatest root of the polynomial. We proceed with the computation of a 5000 
correct digits bound for the positive roots of the polynomial by using a distiller 
whose folding parameter is k = 20 (the respective educated map has order 
m = 2 21 = 2 097 152). We note that this time the graphics are quickly produced 
since we use ListPlot instead of Plot and therefore the respective mesh has 
now much less points than in the usage of Plot. 












Figure 3: Increasing prec to prec = 5000 and for k = 10-folding, the roots on 
[0,1] of the 500-degree Chebyshev’s polynomial are detected. 


4.2 The filtering stage 

Efficient Mathematica commands for manipulating lists, such as Cases, Parti¬ 
tion and Union, are particularly useful in the filtering stage of the distiller 
algorithm. 

Assume that all the numeric points in C have been assigned to a list named data. 
The first step in the filtration deals with the choice of points ( x, y) sufficiently 
close to the bisector line, that is to the line y = x. We test the condition say, 
\y — x\ 2 < 0.1, on the data list, and assign the captured points to a sublist 
named datal, as follows 

datal = Cases [data, x, y; Abs[j/ — x] 2 < 0.1 ] ; (* distille close to 

bissector *) 

The points in datal belonging to a certain platform (that is to a ‘horizontal’ 
segment crossing the bisector line) are good candidates for an approximation of 
the root. The next sublist, named data2, keeps the interesting points. These 
points are chosen in order to satisfy the condition (yi — X\)(y 2 — X 2 ) < 0, which 
assures that at least a machine fixed point exists in the interval denoted by 
[xi,x 2 \: 

data2= 

Cases [Partition[datal ,2,1] ,{{xl_,yl_}>{x2_,y2_} }/; (yl-xl)*(y2-x2)<0 ] ; 

For a given tolerance, say tol = 10 -prec , we are interested in filtering the points 
in the list data2 whose second component y differ from an amount greater than 
tol. Using the command Union, we project the platform ignoring machine- 
duplicate-numbers, obtaining a sublist named union, 

union = Union[Map[Last, Flatten[data2, 1] ]]; 

The final step in the filtering process checks for the accuracy of the former 
captured points. First, one filters the values y in the list union for which f(y) 
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is not greater than say lO -5000 . The result is assigned to a sublist called finalA. 
Second, one filters the values y in the list final A, whose estimated absolute 
error is less than a tolerance, say tol = lO~ 5000 . An error bound for each 
machine root is also computed. The respective code follows. 

finalA = Cases[union, y_ /; f [y] < lO -5000 ] ; 

mapf = Map[{#, g[#, prec] - #} &, finalA]; (* roots and error *) 
tol = 10 —100 ; 

final = Cases[mapf ,{x_, error.}/ ; Abs[error] < tol]; (* error bound 

*) 

Assembling the above filtering stage and the one given at paragraph |2.2[ a 
general function is easily obtainable for the whole distiller algorithm. 

4.3 A 5000-correct digits approximation for the bound of 
the positive roots 

We now apply our distiller to compute a bound for the positive roots of the 
500-degree Chebyshev polynomial in [0.99,1.0], with an error not exceeding 
lO -5000 . Since a bound for the positive roots of the polynomial is required, 
only the greatest computed root in the interval will be displayed as well as its 
estimated error. 

In Figure [4] the ListPlot of the map g is shown, where here g is the Newton’s 
educated method (k = 0), the precision is prec = 5000, and the mesh size 
h = 0.00025. After filtration an empty list is obtained, meaning that this map 
is useless under the previously described filtering criteria. So a more powerful 
‘magnet’ should be used, that is, one needs to increase the order of convergence 
of g by taking a greater value of the folding parameter k. 

_ h=0.00025 40 points _ 


Figure 4: After filtration no points have been captured using the Newton’s seed. 


Increasing to k = 10, the respective 2048-order map g (Figure 5 left) enable us 
to filter relevant points (see Figure [5] right) from which high precision roots can 
be distilled. 
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h=0.00025 29 points 


h=0.00025 41 points 



Figure 5: For a 10-fold map the points in the list data2 of the respective ‘ma¬ 
chine’ step function are shown on the right. 


For h = 0.0025 and k = 20 (Figure 6), the complete filtration process leads 
to 20 fixed points, which are the machine roots of the 500-degree Chebyshev 
polynomial, in the interval [0.99,1.00], for the considered distiller. 


h=0.00025 29 points 


h=0.00025 41 points 1 
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Figure 6: Points captured by a 20-fold g map. 


Denoting by a the last computed root, a 5000 correct digits bound is obtained. 
Respectively the first 100 and the last 100 digits of a are displayed below, as 
well as its estimated error. 


a=0. 

9 999 950 652 018 581 661 118 448 174 487 001 319 149 010 419 592 245 024 005 422 964 693 760 921 636- 
454 0 90 733 615 299 842 527 331. . . 

649 619 039 470 741 295 485 455 419 330 860 833 506 167 181 106 321 342 683 094 894 900 371 436 516- 
812 455 150 525 286 304 558 955 
Estimated error=0. x 10~ 5000 


Note that a Mathematica instruction such as 

x /. Solve[{f[x] == 0,0.99 < x < 1.0}, a;]} 
does not produce an answer within an acceptable CPU running time. 
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Of course, the classical formula giving the zeros of a d-degree Chebyshev poly¬ 
nomial, a.j = cos ^ ^ 2 ’ f° r -i = ■■■(d — 1), can be used in order to confirm 

the above computed value of a. 
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